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Dynamic density-matrix renormalization provides valuable numerical information on dynamic 
correlations by computing convolutions of the corresponding spectral densities. Here we discuss and 
illustrate how and to which extent such data can be deconvolved to retrieve the wanted spectral 
densities. We advocate a nonlinear deconvolution scheme which minimizes the bias in the ansatz for 
the spectral density. The procedure is illustrated for the line shape and width of the Kondo peak 
(low energy feature) and for the line shape of the Hubbard satellites (high energy feature) of the 
single impurity Anderson model. It is found that the Hubbard satellites are strongly asymmetric. 

PACS numbers: 71.55.Ak, 71.27+a, 78.67.Hc, and 78.20. Bh 



I. INTRODUCTION 

The computation of spectral properties is a central is- 
sue in theoretical physics. Many spectroscopic probes 
provide experimental information about the investigated 
systems. In order to understand the meaning of such 
data it is indispensable to be able to compute the cor- 
responding quantities theoretically. This task is particu- 
larly demanding if the system under study is character- 
ized by strong correlations. Then standard approaches 
like diagrammatic perturbation theory have difficulties 
to provide quantitative results. 

An archetypal class of strongly correlated systems are 
impurity models where a small subsystem, the impurity, 
is coupled to a bath of degrees of freedom. The discrete 
levels of the impurity are broadened due to the interac- 
tion with the bath. The most fundamental fcrmionic rep- 
resentative of this class of models is the single impurity 
Anderson model (SIAM) where the impurity is char- 
acterized by a single fermion level which can be empty, 
or singly occupied with spin up or down, or doubly oc- 
cupied. The energy of the doubly occupied state is in- 
creased by the interaction energy U. The bath is a bath of 
non-interacting fermions of either spin direction to which 
the impurity is coupled by hybridization. It is the aim of 
this article to present an algorithm to obtain the dynam- 
ics of a SIAM with high resolution on all energy scales. 

The SIAM describes a plethora of physical problems. 
Historically it was used for diluted magnetic impurities 
in metals, see e.g. Ref. 1. But it also describes the elec- 
tronic transmission through quantum dots, see e.g. Ref. 
2. The smallness of the quantum dot implies a small ca- 
pacitance, hence a large charging energy which represents 
the interaction energy U. The bath is given by the ex- 
ternal leads. The dynamic mean- field theory (DMFT) • 
represents another broad and very active field where the 
SIAM occurs. In this approach, as in all mean-field 
approaches, the lattice problem of strongly interacting 
fermions is mapped onto an effective single-site problem, 
namely a SIAM. This SIAM is linked to the original lat- 
tice problem by a self-consistency condition. The clue is 
that the mean-field, the Green function of the bath, is a 
dynamic quantity depending on frequency. 



The above examples illustrate that it is very impor- 
tant to be able to compute the dynamics of a SIAM in a 
reliable fashion. There arc several numerical approaches 
which can be applied. Among the most prominent ones 
arc quantum Monte Carlo (QMC) 5 and the numerical 
renormalization group (NRG)' 1 ''. Both approaches are 
powerful but do not have a high resolution away from 
the Fermi level. For QMC this is so since the dynam- 
ics is computed in imaginary time and the analytic con- 
tinuation to real frequencies represents an ill-conditioned 
problem. Moreover, care must be taken to reach low tem- 
peratures. The NRG can be used directly at zero temper- 
ature. But it is set up to focus on the limit ui — > 0. The 
energy levels kept are broadened by a broadening which 
is proportional to the frequency which implies that fea- 
tures at higher energies tend to be smeared out 8 . 

We investigate here a third complementary numerical 
approach given by the dynamic density-matrix renor- 
malization (D-DMRG) 8 ' 9 ' 10 ' 11 ' 12 ' 13 ' 14 . Note that there 
is no complete consensus on the nomenclature. Jeckel- 
mann uses the term 'dynamic density-matrix renormal- 
ization group' only for the approach using a variational 
principle 1 while the approach used here is called cor- 
rection vector density-matrix renormalization. We con- 
sider, however, 'dynamic density-matrix renormalization 
group' to be the general term for all algorithms com- 
puting dynamic correlations by DMRG in the frequency 
domain 9 ' 10 ' 11 ' 13 as opposed to DMRG approaches com- 
puting the dynamics in real time 15 ' 16 which efficiently 
implement ideas of Vidal 1 r > 18 . The approach used in this 
article determines the dynamics at zero temperature by 
computing the expectation values in the local propagator. 
This can be realized by targeting not only at the ground 
state and the excited state, but also at the resolvent ap- 
plied to the excited state. This additional targeted state 
is called the correction vector. 

The main limitation of the D-DMRG is that one can- 
not obtain data for purely real frequencies but only for 
frequencies with a certain imaginary part. Hence the 
extraction of the behavior at purely real frequencies is 
one of the main problems to be solved in using the D- 
DMRG 8 ' 14 ' 19 . It is the main aim of the present paper 
to discuss and to compare various algorithms to achieve 
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this extraction. In particular, we will present a nonlinear 
approach from the family of maximum entropy methods. 
This approach provides a continuous, positive ansatz for 
the wanted spectral density with the least bias (LB). 

In the following section II we will present the model 
and we will discuss the observable of interest. Next in 
Sect. Ill, we will discuss various extraction schemes, lin- 
ear and non- linear ones. The features of these schemes 
will be illustrated by some toy spectral densities for which 
the broadened and the unbroadened data is analytically 
available. The results by the LB algorithm for the Kondo 
peak and the Hubbard satellites will be presented in 
Sects. IV and V, respectively. Finally, a short summary 
will be given in Sect. VI. 



II. MODEL AND OBSERVABLE 

We want to illustrate the above general remarks by 
D-DMRG calculations and the subsequent data analysis. 
As motivated in the introduction the single-impurity An- 
derson model (SIAM) at half-filling is a very good and 
interesting testing ground. We will focus on 
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with arbitrary symmetric density of states (DOS) po{u>) 
of the free (U = 0) one-particle Green function Gq(lo) 
of the d-electron. The (i-electron represents the impu- 
rity which is correlated due to the interaction U > 0. 
The bath is represented by the coefficients 7„ > in 
(I). They are the coefficients of the continued fraction 
of the hybridization function T(u>) (cf. Ref. S). Any 
hybridization function with symmetric imaginary part 
Pq(lu) := — 7r _1 ImGo(w + i0+) can be represented by 
an appropriate choice of the j n . Hence the representa- 
tion of the bath as semi-infinite chain does not restrict 
the generality of the model. 

The dynamics wc wish to compute is the dynamics of 
the fermionic single-particle propagator of the d-electron. 
Aiming at the properties at T = it reads 
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Here the ground state is denoted by |0) and its energy by 
Eq . Since we focus at a spin-disordered solution the prop- 
agator has no dependence on the spin index a. Hence, 
it is not denoted as argument of G. The frequencies ui 
and 77 are real. The standard retarded Green function is 
obtained for 77 — ^ 0+ 



Gr(w) = lim G(u> + irf) 
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The quantity we are looking for is the spectral density 
p(u>) := — 7r _1 Im Gr(w). If necessary the real part can 
be obtained from the Kramers-Kronig relation. 

The D-DMRG provides data points at given values of 
lu = £i for finite values of rji 
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where we used the Hilbert representation in the second 
equation. No data can be obtained directly at 77 = since 
the inversion of the Hamiltonian is singular and cannot 
be achieved numerically in a stable way. Henceforth, we 
will call data at finite values of r\ raw data. One way 
to extract the physically relevant data on the real axis is 
to look at a sequence of decreasing values of 77 in order 
to extrapolate the result ' ' to 77 = 0. This approach, 
however, is time-consuming and requires many resources, 
in particular, if one is interested in the whole spectral 
density. So the line followed in the present article is to 
use the raw data as input of a generalized scheme to 
extract the information on the spectral density p(to). 



III. EXTRACTION SCHEMES 

We consider two classes of extraction schemes, linear 
ones and non-linear ones. Linearity means that there is 
a linear relation between the raw data and the wanted 
spectral density. 

A. Linear Extraction Schemes 

1. Deconvolution 

If the unavoidable imaginary part 77 is constant for all 
frequencies Eq. (4) becomes 

1 T • \ 1 f°° rip(uj)du] 

gi = — ImG&+«7 =- / ^ v ' 2 (5) 
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so that the raw data is the convolution of the true spectral 
density p(u>) with the Lorentzian 
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of width 77. Hence the necessary step for retrieving p{io) 
is a deconvolution. It can be achieved in various ways. 
One standard way is to deconvolve the raw data. This 
is done in the time domain reached by Fourier transform 
because the convolution in Eq. (5) becomes a product in 
the time domain 



exp(-J7|r|)p(r) , 



(7) 



where we use p(r) for the Fourier transform of p(to) and 
p(i")raw for the Fourier transform of the raw data. The 
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raw data {gi} is obtained in the first place as a discrete 
set. In order to obtain a quasi continuous distribution 
we interpolate the discrete set {gi} by splines"' : which 
leads to /?(w) raw . The Fourier transforms are most effi- 
ciently done by Fast Fourier algorithms 20 . The actual 
deconvolution is done by dividing by exp(— t)\t\) which 
inverts Eq. (7). Then one transforms back to the fre- 
quency domain. So this procedure is very efficient and 
straightforward* . 

The restriction to be kept in mind is that splining and 
deconvolution cannot create information where no infor- 
mation was present before. If the input data is not pre- 
cise enough or if the spline does not approximate the true 
continuous function p(w) raw well enough the deconvolu- 
tion will fail to produce reasonable results. In practice 
this is seen in unreasonable values of p(r) after the di- 
vision by exp(— 77 1 | ) because this division amplifies any 
inaccuracy for large values of |r|. This problem is cir- 
cumvented by a suitable low-pass filter p T o,At{t) which 
suppresses inaccurate values 

p(T) raw exp(77|r|) -> p(r) raw exp(r?|T|)p TOi Ar(-r) (8) 

at large values of |r| beyond To on the scale At. Of 
course, this implies that only a certain resolution in fre- 
quency can be achieved. In view of the inevitable inac- 
curacies of any numerical calculation one has to accept 
such a bound to the enhancement of the resolution. Nev- 
ertheless, the deconvolution enhances the resolution con- 
siderably and the final curve obtained is continuous for 
all practical purposes since the interpolation allows to 
make the grid as fine as needed. 

The procedure is illustrated in Fig. I where we display 
a power law singularity s(uj) = uj~ a O(uj) exactly and 
convolved by a Lorentzian L 7 (to) of width 7 = 0.05 
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sin (na) (lo 2 + 7 2 ) Q / 2 
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This line shape is well known in photoclectron 
spectroscopy 21 , named the Doniach-Sunjic line shape 22 . 
Note that we changed the normalization so that (9) rep- 
resents exactly the convolution of ui~ a <d(uj) with L 1 (u>). 
The raw data = s(^;7 = 77) that we use for these 
curves is obtained analytically at ^ = i * 0.05 in the 
interval & € [—3,5], see the curve with the circle sym- 
bols. We do not use real D-DMRG here since we want 
first to illustrate the extraction schemes under ideal cir- 
cumstances. The effect of inaccuracies will be discussed 
below. 

In judging the effect of the deconvolution one must 
keep in mind that the reconstruction of a diverging sin- 
gular line is the worst case for any algorithm. We have 
chosen the Doniach-Sunjic line shape for illustration in 
order to highlight the differences in the various schemes. 
Below (Fig. 3) we will present results also for a smoother 
curve to show that such a curve can be reconstructed in 
a quantitatively reliable way. 
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FIG. 1: Doniach-Sunjic line shape (Eq. (9)) for a = 1/4 
with (thin black solid line with circles representing the raw 
data at = 0.05i) and without broadening (thick black solid 
line) 7 = j) — 0.05. The results of various schemes to re- 
trieve the unbroadened line are shown: by deconvolution via 
fast Fourier transform (FFT), by matrix inversion assuming 
spikes (MI(S)) or histograms (MI(H)) or a piecewise linear 
continuous function (MI(C)), by the non-linear least-bias al- 
gorithm (LB). All schemes are described in the main text. 



The dashed-dotted line is the result of the above de- 
scribed deconvolution using the low-pass filter 
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exp tan y ^kJ" J } otherwise 

(10) 

which was also used to achieve the deconvolutions shown 
in Fig. 2 of Ref. 8. The parameter To determines where 
the low-pass cutoff is done; the parameter At determines 
on which time-scale the cutoff function switches from 1 
to 0. We used To « 10.7 and At f=a 0.763. In practice, it 
turns out that it is fairly obvious in the T-domain which 
values one has to choose for the low-pass filter. The data 
for too large values of |t| scatter very much. 



2. Matrix Inversion 

A robust alternative to the deconvolution by Fourier 
transform is the explicit matrix inversion of the convolu- 
tion procedure. This procedure shares the linearity with 
the Fourier transform and it may also lead to negative 
spectral weight close to abrupt changes of p(ll>), for in- 
stance at singularities. An advantage over the Fourier de- 
convolution is that one may also consider variable widths 
r/i. In principle, this allows to adapt the grid {^} to the 
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FIG. 2: Illustration of various assumptions on the behavior 
of the spectral density p(u). (a) Spikes: p{ui) consists of a set 
of (^-functions with given positions, but unknown weights, (b) 
Histogram: p(ui) is piecewise constant with given positions of 
the jumps, but unknown heights pi. (c) Continuous: p(ui) is 
piecewise linear and continuous with given positions of the 
cusps, but unknown values p(u>t). 



expected behavior of p{ui). A denser grid can be taken 
where the spectral density varies more rapidly. 

The raw data {gi} provided by D-DMRG for a set of 
are taken as the components of the vector g. Then 
the process of convolution can be described by a linear 
mapping M of a set of linear parameters {k} character- 
izing p(u>) onto the raw data {gi}. Let us take the {li} 
also as components of a vector 1. Then the convolution 
reads 

g = Ml. (11) 

Clearly, the deconvolution implies the inversion of this 
equation. To this end, we have to choose the matrix M to 
be square which means that there must be as many data 
values gi as there are parameters li to be determined. It 
is obvious that this analysis is linear. The precise form 
of the matrix M depends on further assumptions. Three 



generic scenarios are studied; they are illustrated in Fig. 
2. 

a. Spikes Assuming that p{u>) is given as set of 5- 
functions (cf. Fig. 2a) with weights {wi} at given fre- 
quencies {uJi}. Then the weights constitute the linear 
parameters defining the spectral density. The matrix M 
is derived from Eq. (4); it has simple matrix elements 
given by Lorcntzians 

M nii = Lr^dn - Ui) . (12) 

This approach has been proposed and used by Jeckel- 
mann and coworkers 11,10 . To have a square matrix prob- 
lem the number of spikes has to equal the number of gi 
determined by D-DMRG. 

The rendering of the results is subtle. Since one cannot 
plot (5-functions they have to be broadened which adds 
another free parameter. Which broadening one has to 
choose is not a priori clear. For equidistant grids 
the distance between two consecutive peaks is a natural 
choice 14,11 ', cf. Fig. 1. 

b. Histogram One assumes that p(w) is piecewise 
constant p(ui) = pi for LOi < u> < cjj+i (cf. Fig. 2b) where 
the frequencies u>i are given beforehand. The set of linear 
parameters defining the spectral density are the values 
li = pi. Again, Eq. (11) has to be solved. The matrix 
elements of M are found from the integration of the right 
hand side of Eq. (4) 

M„s = tt" 1 arctan [(cj - £ n )/r/ n ] . (13) 

The number of values pi (not the number of £j,) must be 
equal to the number of & to ensure that M is square. The 
rendering is straightforward since the values pi represent 
densities which can be plotted directly, cf. the histograms 
in Fig. 1. 

c. Continuous One assumes that p(u>) is continuous 
and piecewise linear (cf. Fig. 2c) with p(u)j) = Pj,j £ 
{1, 2, 3, . . . ,p} at given frequencies u>i. The values pi rep- 
resent the linear parameters li — pi which are determined 
by Eq. (11). The matrix elements result from the inte- 
gration in Eq. (4). They are given by M n ,i — d Pi A n (£ n ) 
with 

p-i 

MO = E {^i Hvl + (£ - ^) 2 ]/(2tt)+ 

3 =i (14a) 

(£aj + b j) arctan[(w - 0/Vn]M 
a j = {Pj+i - - (14b) 

h j = (uj+iPj - Vjpj+i)/ (Wj+i - cjj) . (14c) 

The number p must be equal to the number of raw data 
points. An example is depicted by the piecewise linear 
dashed curve in Fig. 1. 

In all schemes, the numbers must be chosen such that 
M is a square matrix. This is a necessary but not a suffi- 
cient condition for the existence of a unique solution 1 in 
Eq. (11). In practice, we did not encounter problems in 
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the inversion of Eq. (11) as long as the raw data points 
were distributed rather evenly along the real axis. Only 
if there are data points accumulating in certain regions, 
for instance several data points at the same frequency, 
the inversion can be problematic. Loosely speaking this 
may occur since the raw data is slightly contradictory 
due to numerical inaccuracies. The broadening by the 
{77™} reduces the differences between the spectral densi- 
ties. Hence small deviations in the raw data have large 
effects on the extracted spectral densities. 

All linear extraction schemes do not guarantee that 
the extracted spectral density is non- negative, see Fig. 1. 
Whether this must be considered a serious drawback de- 
pends on the extent to which negative values occur and 
on the context in which the result is used. If the spectral 
density is the final result small regions of overshooting 
are unproblcmatic. If, however, the overshooting is con- 
siderable and if the spectral density shall be used in a 
subsequent step, e.g. in the self-consistency of a DMFT 
calculation, then negative values pose a severe problem. 
An important example is the determination of the coef- 
ficients of the continued fraction of the spectral density. 
This determination is only possible if the spectral density 
is really non-negative. 

Another problem is the smoothness of the extracted 
density. Spurious discontinuities like the ones assumed in 
the spike ansatz or in the histogram ansatz (see Figs. 1 
and 2) imply singularities in the real part of the propaga- 
tor which is determined by the Kramers-Kronig relation. 
This in turn leads to unwanted features like slowly decay- 
ing oscillations in coefficients of the continued fraction. 
Of course, various schemes can be used to interpolate 
the discrete data provided by the matrix inversion ap- 
proaches. But the interpolation represents an additional 
approximation which can be difficult to control. 



B. Nonlinear Extraction Schemes 

1. Basic Algorithm 

In view of the drawbacks of the linear extraction 
schemes it is worthwhile to think about alternatives. The 
objective is to devise an ansatz for a continuous, non- 
negative spectral density p(w) which is consistent with 
the numerically determined values of the raw data {gi} 
in Eq. (4). The ideal ansatz is completely unbiased. That 
means it does not use any information other than the one 
provided by the raw data. The information content of a 
density p(u>) is measured up to a constant by its negative 
entropy 
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(15) 



The least biased ansatz is the one with the least informa- 
tion content which is still compatible with the raw data. 
Hence we have to look for the density p(u>) which mini- 
mizes —S (maximizes S) under the conditions (4) given 



by the raw data {gi} and by the known normalization 

/oo 
p(uj)duj . (16) 
-00 

To find this least biased ansatz (LB) is a straightfor- 
ward task. Using the Lagrange multipliers {Xi} for the 
p conditions set by the raw data {gi} and the Lagrange 
multiplier p for the normalization (16) the least biased 
ansatz is characterized by SS = 0, i.e. 

p 

= -l-]n(p(w)) + 53AiA»(w-&)+A* ■ ( 17 ) 

i=l 

This equation implies that the LB ansatz reads 



p(u>) = exp 
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(18) 



where we defined p = p — 1. The Lagrange multipliers 
are determined by the nonlinear equations (16) and (4). 
They can be determined by any standard algorithm for 
a set of nonlinear equations. Via the ansatz (18), the 
p + 1 Lagrange multipliers determine the most unbiased 
spectral density p(uj) which is still compatible with the 
numerically measured information on p(u>). 

The LB ansatz (18) is positive and continuous. Hence 
it avoids two major drawbacks of the linear extraction 
schemes. In spite of its continuity the LB ansatz is gov- 
erned by a restricted number of parameters. Despite 
its positivity, the LB ansatz is able to reproduce rather 
abrupt changes in the spectral density, see Fig. 1. If 
arbitrarily accurate data at a fixed value of 77 were avail- 
able on an arbitrarily dense grid, the LB algorithm were 
able to provide the correct result with arbitrary resolu- 
tion (see also the discussion of Fig. 4 below). But the 
accuracy of the raw data required to achieve a certain 
resolution in the deconvolved result grows exponentially. 
So in practice this route cannot be followed very far and 
the broadening r\ sets the scale for the achievable resolu- 
tion. 

The least bias approach belongs to the class of max- 
imum entropy methods (MaxEnt)-". The main differ- 
ence to standard MaxEnt is that we do not use a x~ 
functional in addition to the entropy function (15). The 
X-functionals are bilinear in the density. They are intro- 
duced to account for possible deviations of the from 
their true values. Such deviations occur for instance in 
quantum Monte Carlo calculations due to the inevitable 
statistical error. The D-DMRG data is free from sta- 
tistical errors. Hence we can use the entropy functional 
alone as described above. The correction of numerical 
inaccuracies is discussed in more detail below. 

We emphasize a major difference between the extrac- 
tion of the spectral density p(u>) from D-DMRG data and 
from QMC. In the former case the task to be solved is to 
remove a small imaginary part u> + if] — > u>. This means 
that in the D-DMRG the raw data is situated slightly 
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above the real axis and has to be continued down to it. 
In the latter case the task to be solved is to continue 
QMC data from the imaginary axis (no real part) to the 
real axis iC, — * uj. So the challenge in the QMC data 
analysis is much greater than in the D-DMRG data anal- 
ysis. This explains why the D-DMRG approach is suited 
to investigate sharp features also at high energies s while 
this is not a straightforward task by QMC. 

The solution of the nonlinear equations (16) and (4) 
is done by a standard algorithm. There is no mathe- 
matical argument to show that there is only one unique 
solution. For instance, we found that the normalization 
of the spectral density need not be ensured by the pa- 
rameter p. This parameter can be fixed to almost any 
value. The remaining Lagrange parameters suffice to de- 
termine good approximations to p{uj) which fulfill the 
normalization condition (16) well. This implies that it- 
erative numerical solutions have difficulties to fix p in- 
dependently. The root-finder algorithms run much more 
stable if the normalization is not included in the set of 
equations. Yet the resulting densities are normalized if 
the raw data provides a reasonable scan of the spectral 
density, i.e. if there is raw data at all frequencies uj where 
the density is non-negligible. 

For other sum rules, for instance the second moment 
(<jj 2 ) , the same conclusion holds. If the raw data scans all 
relevant frequencies the sum rules do not provide useful 
additional information. But if raw data is only available 
for restricted frequency intervals, the sum rules help to 
improve the LB ansatz. 

Furthermore, the iterative numerical solutions depend 
sometimes on the initial values. But the resulting p(uj) 
are in general (almost) identical. This remains true if the 
iterative algorithm does not find a true solution of the 
set of non-linear equations but only a set of parameters 
which makes the deviations 

Agi := gi / T rj- — 2 ( 19 ) 

small but fails to make them zero. In summary, the iter- 
ative determination of the Lagrange multipliers does not 
represent a major problem. 

In Fig. 1 we display a comparison of all data extraction 
schemes introduced so far. The line shape to be found is 
a singular power law divergence. This line shape consti- 
tutes an unsolvable task since a divergence cannot be re- 
produced by the algorithms discussed. But this example 
illustrates well to which extent the algorithms manage 
to render the true distribution of spectral weight. All 
the linear schemes lead to regions of negative spectral 
densities which is a severe drawback. The divergence is 
approximated by a broad peak located at 0.05 away from 
the position of the divergence. There are some spurious 
oscillations in the approximated spectral density. 

The LB scheme avoids negative spectral weight by con- 
struction. The divergence is approximated by a sharper 
peak at about 0.02 away from the position of the diver- 
gence. The algorithm implies spurious oscillations in the 
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FIG. 3: Exact (thick black solid line) and broadened (by 
rj = 0.05, thin black solid line with circles representing the 
raw data at £j = 0.05i) lines derived from the Doniach-Sunjic 
line shape in Eq. (9), see also main text. The results of var- 
ious schemes to retrieve the unbroadened line are shown: by 
deconvolution via fast Fourier transform (FFT), by matrix 
inversion assuming spikes (MI(S)) or histograms (MI(H)) or 
a piecewise linear continuous function (MI(C)), by the non- 
linear least- bias algorithm (LB). 



approximated spectral density. So we conclude that the 
LB analysis represents a very efficient reconstruction of 
the unbroadened data even though the spurious oscilla- 
tions can be a nuisance. 

On the other hand it is, of course, possible to improve 
the analysis. An additional data point at w w 0.02 will 
certainly help all algorithms to reproduce the unbroad- 
ened density more faithfully. The same is trivially true 
for a denser mesh of raw data points. The former so- 
lution, however, requires either to intervene manually in 
the data analysis or to know beforehand where the peaks 
will be located. The latter requires much more numerical 
effort on the D-DMRG level so that this is not the most 
efficient approach. 

In Fig. 3 we present the analysis of a smooth exact 
curve, namely the line shape in (9) at 7 = 0.01. The 
curve broadened by 77 = 0.05 is the one at 7 = 0.06 since 
the widths of Lorentzians is additive under convolution, 
in contrast to the root-mean-square of narrower distri- 
bution functions. Clearly, the extraction schemes do a 
better job for this non-singular case. The regions of neg- 
ative spectral density in the linear schemes has shrunk. 
The LB scheme manages to reproduce the true density 
almost perfectly. The spurious oscillations arc negligible. 
If we had chosen 7 = 0.02 for the exact curve the LB 
density would be hardly distinguishable from the exact 
curve. 
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FIG. 4: Exact (thick black solid line) and broadened (by 
rj = 0.05, thin black solid line with circles representing the 
raw data at £j = 0.02i) lines derived from the Doniach-Sunjic 
line shape in Eq. (9). The dashed line depicts the density 
extracted from exact raw data by the LB scheme (the curve 
lies on top of the thick black line except at the maximum). 
The dotted oscillations are the LB result from data contam- 
inated by a random error, see main text. The thin solid line 
represents the density derived from the contaminated data by 
the robust LB scheme with A = 1 using Eqs. (21,22). 



2. Robustness towards Inaccuracies 

So far we analyzed ideal raw data, i.e. no errors were 
considered. Statistical errors do not occur in the D- 
DMRG approach but systematic errors occur. There are 
two main sources for such errors. The first is the inaccu- 
racy of the algorithm due to the truncation of the Hilbcrt 
space. This is an unavoidable error; but it can be con- 
trolled by comparing the results for different numbers of 
states kept in the density-matrix renormalization. We 
perform our calculations for 128 and for 256 states us- 
ing the representation of spinful fermions as two kinds of 
spins via a double Jordan- Wigncr transformation '. The 
relative truncation error in the spectral densities is esti- 
mated to be of the order of 10 -5 to 10~ 3 depending on 
the frequency where it is computed. For low values of the 
frequency (\uj\ ;$ D in our model) the lower error applies; 
for frequencies beyond D the larger value applies. 

The second important source of inaccuracy are finite- 
size effects. In principle, we wish to compute the spectral 
density for the infinite system. But this is not feasible 
numerically. So the system - the infinite chain - is ap- 
proximated by finite chains of L = 120 to 400 sites. In 
a rigorous sense, the spectral density of the finite system 
consists of (*>- functions, i.e. it is not a continuous func- 
tion. But we do not intend to resolve all the fine details 
of the finite system. Rather we interpret our numerical 



raw data obtained for the finite chain as an approximate 
description of the infinite system. The deviation of the 
raw data for the finite chain from the desired raw data of 
the infinite chain is considered the source of a systematic 
error, the finite-size effect. There are two conceivable 
ways to deal with this error. One way is to perform an 
extrapolation in system size L for each raw data point 
gi before the deconvolution is done. The other way is to 
deconvolve the raw data for various chain lengths and to 
check whether the results still depend on the length L . 

In our work, we have chosen the second approach be- 
cause the first is hampered by an unsystematic behavior 
of gi on L. Depending on details a particular <5-pcak of 
the finite systems contributes more or less to the gi under 
study. This makes a controlled extrapolation for all {gi} 
difficult, if not impossible. 

In the second approach, care must be taken that the 
length L is so large that the raw data is sufficiently close 
to the raw data of the infinite system. In practice, this 
puts a restriction on rj and L, see Sect. IV and Ref. 13. 
Of course, the use of finite chains restricts the extent to 
which we may extract information on the exact infinite 
system from the broadened data obtained for the finite 
system. 

Another sort of errors are rounding errors. But they 
arc of minor importance compared to the two other 
sources discussed above. 

In Fig. 4 we display the LB analysis of raw data for the 
exact curve at 7 = 0.01. In contrast to the procedure in 
Fig. 3 we use a finer grid of & = i* 0.02 for the raw data 
(curve with circles). The extracted curve represents the 
exact one very well, see the dashed line in Fig. 4. The 
agreement is significantly better than the one reached in 
Fig. 3. This illustrates that sufficiently accurate broad- 
ened data at fixed rj can be used to resolve features of 
widths below r\. 

As explained above, there are in practice restrictions 
to better resolutions due to the systematic errors, namely 
the truncation of the basis and the finite-size effect. To 
examine the effect of systematic errors on the LB decon- 
volution we deliberately contaminated the raw data by 
an error of the order of 10 -3 according to 

9l 9i * (1 + 10~ 3 x) (20) 

where x is a random number between —1 and 1. The 
randomness is just used to mimic a systematic error 
which is uncorrelatcd from frequency to frequency. We 
obtained qualitatively very much the same results for a 
non-random error gi — > gi * (1 + 10 -3 cos(v / 5*))- 

If the data is slightly inaccurate the deconvolution in- 
deed fails as illustrated by the wild oscillations of the 
dotted line. We conclude that the occurrence of strong 
oscillations can be taken as criterion that the used raw 
data is not accurate enough for the LB analysis, i.e. the 
systematic errors are too large. 

The thinner solid line depicts a successful deconvolu- 
tion of the contaminated raw data. It is achieved by a 
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modification of the LB algorithm which makes it a stan- 
dard maximum entropy approach. The negative entropy 
functional (15) is supplemented as shown 

F[p(u>)}: = -S[p(u)} + A X [p(u)] (21) 

/OO 
pHln( (0 M)d W + A^(A ft ) 2 
-oo • 

by a quadratic functional \ which measures the distance 
from the perfect fulfillment of the conditions (4). The 
differences Agi are defined in Eq. (19). The minimization 
of the functional F leads to the same ansatz (18) as before 
except that the parameters Aj are now given by 

A; = 2AA 9l . (22) 

The set of these nonlinear equations (instead of Agi = 
as for the pure LB approach) is used to determine the 
parameters A;. It is obvious that the robust modification 
of the LB ansatz becomes the pure LB ansatz in the limit 
of A — > oo since in this limit the deviations Agi vanish 
for given values of the Lagrange parameters A^. If the 
data is only weakly contaminated by inaccuracies, large 
values of A can be used to extract the spectral densities. 
Our example in Fig. 4 shows fairly strongly perturbed 
data. Still the robust LB can make sense out of them 
and retrieves a good approximation to the underlying 
curve. 

The modified LB ansatz (21) is more robust since it 
can deal with some inaccuracies or inconsistencies of the 
raw data. Imagine that we deal with raw data on a fine 
grid where the distance between the data points is signifi- 
cantly smaller than the Lorentzian width r\i 3> (£;+i — 
Then gi+i and gi may differ only slightly if they are de- 
rived from a smooth continuous density p(to). Any in- 
accuracy spoils this relation and may introduce signif- 
icant spurious oscillations, sec Fig. 4. As we stressed 
already previously the broadening makes different data 
more alike. Hence, the inverse process enhances slight 
differences like the ones between exact and inaccurate 
raw data greatly. The robust LB ansatz (21) helps to 
make the data extraction less sensitive to such effects 
without losing much resolution. Thereby, spurious oscil- 
lations can be suppressed. 

The robust LB scheme opens the possibility to resolve 
features of widths below a given value of the broadening 
r\ by using a finer grid with A£ < rj since slightly inaccu- 
rate data can still be deconvolved. In this way, one may 
avoid the explicit use of small values of r\. We emphasize, 
however, that the data must be accurate enough to con- 
tain the information on the relevant physics. Of course, 
the robust LB approach is no means to extract informa- 
tion which is not given by the raw data. For instance, one 
may stick to short chain lengths only if the underlying 
physical problem does not demand to describe long-range 
spatial fluctuations. 
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FIG. 5: Kondo peaks for the single-impurity Anderson model 
(1) for a hybridization function with semi-elliptic spectral 
density and V = D/2. On increasing interaction the Kondo 
peak becomes rapidly narrower. 



IV. KONDO RESONANCE 

Next we present data for the local spectral density p{u>) 
for the SIAM as introduced in Sect. II. In the present 
section we focus on the behavior at low energies where 
p(u) is dominated by the Kondo peak, see e.g. Rcf. 1. 
This peak can be computed by NRG in a very efficient 
way 1 ' 6 ''. So the idea of the present calculation is not 
to provide novel data, but to gauge the D-DMRG and 
to demonstrate that features at low energies can be re- 
solved. We use the LB scheme to extract the real spectral 
density from the raw data as explained above. This is 
one of the differences to the investigation by Nishimoto 
and Jeckelmann 14 . For details of the numerical realiza- 
tions we refer to previous works 8 ' 14 . Another difference 
concerns the parameter regime. We consider here inter- 
action values U beyond the bare band width W = 2D, 
i.e. U > 2D, while Nishimoto and Jeckelmann look at 
Lorentzian bare spectral function whose band width is 
very large. But for the low energy region this difference 
is only a quantitative one. 

Our results are depicted in Fig. 5. We have chosen the 
parameters such that the hybridization function has a 
semi-elliptic spectral density. The hybridization is taken 
to be V = D/2 so that the spectral density of the d-site 
without interaction is semi-elliptic, too. Clearly visible 
is the rapid narrowing of the Kondo peak on increasing 
interaction. Note that the density at zero energy p(0) is 
pinned to its non-interacting value as required by the sum 
rules 1 ' 2,5 ' 24 ' 25 . This exact result is fulfilled to numerical 
accuracy which ranges from 0.1% for smaller interactions 
to 1% at larger interactions, see Fig. 5. 

The data in Fig. 5 is obtained on various grids in fre- 
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quency, for various chain lengths and for varying widths 
r). All calculations are done with m = 128 sites kept in 
the truncations of the DMRG. The chain lengths L vary 
between 120 and 400 fermionic sites which translates to 
240 (e.g. at U = 1.25D) to 800 (e.g. at U = 4.00D) 
spin sites after the double Jordan- Wigncr transforma- 
tion. The smaller the widths 77 arc chosen the longer the 
chains have to be taken 12 ' 13 . A chain of L fermionic 
sites implies L main peaks distributed over the band 
width W = 2D. Assuming a roughly equidistant dis- 
tribution the distance between two main peaks is W/L 
which should be significantly smaller than 77 



2W/L < 77 



(23) 



in order to ensure that the discrete structure of the finite 
system is sufficiently smeared out. Then the data pro- 
vided can be interpreted reasonably well as data of the 
infinite system. 

Since the calculations are less costly at low energies 
than at higher energies we used mixed raw data coming 
from various chain lengths. The width r\ is varied cor- 
respondingly; we always used A£ = — £j = rj. For 
instance at U = 4D, we used L = 400 with 77 = 0.01Z? 



between u> 
between u> 
between u> 



0.00 and 0.05 
0.06 and 0.18 
0.20 and 2.95 



L = 200 with 77 = 0.02D 
L = 200 with 7/ = 0.05D 
L = 120 with 7/ = 0.10D 



between lu = 3.00 and 4.00. The analysis of the raw data 
was done in all cases by the pure LB ansatz (18) with the 
conditions Agi = 0. There was no need to use the robust 
LB with Aj = 2AAgi. 

The rapidly narrowing peaks in Fig. 5 are characterized 
by the Kondo energy scale Xk- This scale can be read off 
from the spectral densities, for instance as half the width 
at half the maximum, i.e. at TrDp(uj = Tk) = 1. From 
analytic considerations 1 , it is known that the Kondo en- 
ergy scale is exponentially small in the interaction 



T K oc V^D/U cxp (-TiUD/{Wf 



(24) 



This formula applies in the limit of U W = 2D which 
is called the Schrieffer- Wolff limit of the SIAM. It holds 
for U — ► 00. The numerical results for the Kondo scale 
and the analytical prediction (24) agree very well in Fig. 
6. Only for U < W deviations occur as was to be ex- 
pected. The widths at large values of U (U > AD) arc 
found from results for p(ui) where we extracted only the 
behavior at low energies. To determine the HWHM it is 
not necessary to know the whole line shape. It is an asset 
of the LB extraction scheme that it allows also to deter- 
mine only a part of the whole curve. In conclusion, Fig. 6 
demonstrates that the dynamic DMRG is able to repro- 
duce the low energy scale of the SIAM over two orders of 
magnitude. 



V. HUBBARD SATELLITES 

In this section wc show that the D-DMRG combined 
with the powerful LB scheme allows to resolve features 
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FIG. 6: Kondo energy scale (Kondo temperature) for the 
single-impurity Anderson model (1) for a hybridization func- 
tion with semi-elliptic spectral density and V = D/2 in the 
Schrieffer-Wolff limit U > 2D. Symbols: half-width-half- 
maximum (HWHM) read off from p(uj) as found by the LB 
analysis of D-DMRG raw data; line: analytic result in Eq. 24. 
The proportionality factor of the fit is cm = 3.887. 
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FIG. 7: Upper Hubbard satellite for two values of the inter- 
action U obtained from D-DMRG raw data (rj = 0.1D, A£ = 
0.05D, m = 256, L = 120) processed by the LB scheme. At 
larger U a strongly asymmetric line shape occurs. 



at high energies which so far eluded a quantitative de- 
termination. The calculation of the line shapes of energy 
levels at high energies represents a new field of applica- 
tions since previous methods are not suited to perform 
such computations. An exception are features at high en- 
ergies which can be understood and described as shifted 
low energy features 21 '. We emphasize that the progress in 
the manipulation of quantum dots by optical means has 
brought the measurements of high energy features within 
reach, see e.g. Refs. 27 and 28. 

In a previous work s , we have shown that the Hubbard 
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satellites at about u> — ±{7/2 become increasingly sharp 
once the interaction U is larger than the bandwidth W = 
2D. In that work, however, we were not able to resolve 
the sharp Hubbard satellites since the deconvolution by 
Fourier transform was not powerful enough to do so, cf. 
Sect. III. By means of the LB extraction we are now 
able to address the line shape. In Fig. 7 the results for 
a moderate value of U (2.57)) and for a larger value of 
U (3.57)) are shown. While the line at U = 2.57) is still 
fairly symmetric (cf. also the line for U = 2D in Fig. 2 
of Ref. 8) the line at U = 3.57) has a clear asymmetric 
shape. The rise at the low energy side is rather abrupt 
and steep while the decrease at the high energy side is 
much more gentle and slow. The peak is very pronounced 
and the maximum value very high. The small bumps left 
of the peak are artefacts of the data extraction similar 
to the small dip in front of the uprise in Fig. 3. We also 
consider the shoulder on the right side of the peak in Fig. 
7 to be an artefact of the data analysis in analogy to the 
spurious oscillations in Fig. 1. 

At present, we do not have an analytical description 
of the line shape in Fig. 7. Qualitatively, however, the 
following description holds. The doubly occupied level 
has the energy U/2 in the atomic limit for V = 0. This 
level is shifted to higher energies by a finite hybridization 
V > 8 ' 29 and it is broadened*. The shift oc V 2 /U is due 
to level repulsion between the doubly occupied level and 
modes where the doubly occupied level excites additional 
particle-hole pairs. A rough, averaged description of the 
broadening is provided by Fermi's golden rule implying 
that the width is proportional to V 4 /U 2 . 

In view of the asymmetric line shape and of the very 
narrow peak this averaged picture can be improved. The 
narrow peak itself has a certain intrinsic width which is 
small. Note that the general understanding of the SIAM 
as a local Fermi liquid implies that there are no singular- 
ities in the propagator. The spectral density is a smooth 
function, though it may display very narrow features. 
This is in contrast to the X-ray edge problem where the 
added fermion is completely immobile (no recoil). 

The intrinsic width of the narrow peak results from the 
decay of the doubly occupied level into particle-hole pairs 
of the local Fermi liquid. According to Fermi's golden 
rule such a decay can take place only if the energy of the 
initial and the final state is equal which means that many 
particle-hole pairs have to be excited in order to make 
up for the relatively large energy m U/2 of the double 
occupancy. Very many particle-hole pairs are needed in 
view of the reduced effective band width of the order of 
the Kondo energy scale 7k . This argument accounts for a 
finite, though very large, life time of the doubly occupied 
level as seen in Fig. 7. 

The high energy tail of the line in Fig. 7 can be under- 
stood as an effect of particle-hole pairs which are gener- 
ated by the double occupancy. This is not the same as 
the decay of the double occupancy described above. In 
the decay the double occupancy disappears in the pro- 
cess while it remains on generating additional particle- 



hole pairs. This explains why the contribution of this 
process is found at the high energy side of the peak in 
Fig. 7. The additional particle-hole pairs require an ad- 
ditional amount of energy to be created. The physics of 
this process is similar to the physics of the X-ray edge 
problem where the change of a local potential induces 
infinitely many particle-hole pairs thus leading to slowly 
decreasing tails in the spectra, see for instance Ref. 30 
and references therein. The main difference to the X-ray 
edge problem is that there is no true singularity here so 
that the line shape is smoother. 

At present, the accuracy of the data in the high energy 
tails is not sufficient to search for approximate power laws 
which possibly describe the tails of the Hubbard satellites 
in Fig. 7. 



VI. SUMMARY 

In this article, we discussed a variety of schemes to 
extract the spectral density p(u>) from the results of dy- 
namic density-matrix renormalization data. All these 
schemes have the aim to remove the unavoidable broad- 
ening which has to be included in a D-DMRG calcula- 
tion. The linear schemes use either Fourier transform to 
deconvolve the raw data or they implement an explicit 
matrix inversion. These schemes are linear because there 
is a linear relationship between the raw data and the ex- 
tracted spectral density. If the structures to be resolved 
are not too sharp the linear schemes work well. If there 
are sharp structures the linear schemes are prone to lead 
to negative spectral densities which result from spurious 
oscillations. Furthermore, they can resolve the positions 
of sharp peaks only with the accuracy of the grid on 
which the raw data has been computed. 

The nonlinear scheme introduced belongs to the family 
of maximum entropy methods. If the raw data is suffi- 
ciently accurate the least-bias approach works very well. 
It provides a positive and continuous ansatz for the spec- 
tral density with the least possible bias. Even relative 
abrupt changes of the spectral density can be reproduced 
satisfactorily. In the vicinity of singularities spurious os- 
cillations occur. But they do not violate the positivity 
of the ansatz. The least-bias ansatz can be made more 
robust towards small numerical inaccuracies and finite- 
size effects by including besides the entropy functional a 
X-functional in the functional to be minimized. Thereby, 
one can allow for small deviations from the raw data. 

The properties of the above mentioned schemes were 
illustrated by calculations for a singular toy spectral den- 
sity of the Doniach-Sunjic type. D-DMRG calculations 
were carried out and presented for the Kondo energy scale 
of the symmetric single impurity Anderson model in the 
Schrieffer- Wolff limit of U > W where U is the interac- 
tion and W the band width. The Kondo energy scale 
could be retrieved over two orders of magnitude. 

The reproduction of the Kondo energy scale served as a 
gauge for the approach. To our knowledge, numerically 
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exact results like those obtained for the line shape of 
the Hubbard satellites at uj = ±{7/2 beyond the region 
of the bare band ui E [—D,D] are not available in the 
literature. By these line shapes we extended our previous 
analysis of the Hubbard satellites s . Furthermore, these 
line shapes illustrate that the dynamic density-matrix 
renormalization is suited to provide high resolution data 
at high energies, i.e. away from the Fermi level. It is 
to be expected that such information will become more 
and more important as the experimental techniques are 
improving. 
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